Short-range interacting system without additivity 
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It is argued that a certain kind of short-range interacting systems exhibits nonadditivity when 
several timescales are well separated. Under the condition of separated timescales, the system is 
described by the elastic spin model. We numerically find that it is extensive but nonadditive, 
which is directly confirmed by the work measurement and also indicated by ensemble inequivalence. 
Further, we estimate the effective Hamiltonian for the spin variables, and it is clarified that the 
effective interaction is long-ranged. Remarkably, so called Kac's prescription, which is usually 
regarded as a mathematical operation to make the system extensive, naturally holds. 



Let us consider a system consisting of two macroscopic 
subsystems A and B. In a usual macroscopic system, 
the total amount of energy is given by sum of internal 
energies of the two subsystems because the interaction 
energy between A and B is negligible compared to the 
bulk energy. This property is called additivity (the pre- 
cise definition will be given later) . Additivity is regarded 
as a fundamental property of macroscopic systems [1]. 
It ensures concavity or convexity of the thermodynamic 
function. In statistical mechanics, it leads to ensemble 
equivalence, i.e. several statistical ensembles yield iden- 
tical thermodynamic quantities [2]. However, not all the 
macroscopic systems possess additivity. Long-range in- 
teracting systems are representative of nonadditive and 
physically relevant systems [3, 4]. Because of the lack 
of additivity, long-range interacting systems can exhibit 
unfamiliar and peculiar macroscopic properties, e.g. the 
negative specific heat [5, 6], ensemble inequivalence [5], 
macroscopic inhomogeneity [7] , and no thermalization in 
an isolated system [8] . 

Apparently, a short-range interacting system is un- 
likely to be nonadditive since the interaction energy will 
be very small compared to the bulk energy. In this Letter, 
however, it is pointed out that there exists a short-range 
interacting system without additivity. Interestingly, in 
spite of its nonadditivity, the model is extensive; the en- 
ergy is proportional to the system size if we make the 
system large uniformly. 

Now we explain the "elastic spin model" studied in 
this Letter. The model itself has been already known 
in studies on spin-crossover materials [9]. See Ref. [10] 
and references therein for detailed explanation of spin- 
crossover materials. N molecules are aligned on the 
two-dimensional triangular lattice of side L, where N = 
L 2 .[18] Each molecule is composed of a metal ion and 
surrounding ligands. Each molecule i has two different 
stable internal states cr, = ±1, that is, the high-spin (HS) 
state <7j = +1 and the low-spin (LS) state at = —1. Elec- 
tron configurations in a metal ion are different between 
the HS and the LS state and the HS state has a higher 
spin than the LS state. For example, in Fe 11 , the HS state 
has S = 2 and the LS state has 5 = 0. As a result, the HS 
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FIG. 1: (a) A typical shape of the interaction potential, (b) 
Triangular lattice and the label of the molecules. 



state with S — 2 has degeneracies of S z — 2, 1, 0, —1, —2. 
In general, degeneracies of the HS state are denoted by 
a parameter g and thus <7j G {+1, +1, +1, • • • , +1, — 1}. 
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Although Oi is not a genuine spin, we call it a "spin 
variable" and M = G i the "magnetization" . The 
"magnetization density" is denoted by m = M/N. 

The intermolecular interaction is given by some short- 
range potential Vy(r) such as Fig. 1(a), which decays 
faster than 1 jr 2 in a long distance. Because of the short- 
range nature of the interaction, it is sufficient to consider 
only the nearest-neighbor interactions. The important 
point is that the equilibrium distance of Vij{r) is given 
by Rij = R((7i) + R((Tj), which depends on the molec- 
ular internal states. The values of and R(— 1) 
represent the molecular radius at state <7j = +1 and — 1, 
respectively. This size difference is actually observed in 
experiments [11]. It has been pointed out that this size 
difference plays an important role for spin-crossover tran- 
sitions [12]. 

When the potential depth Vb in Fig. 1(a) is much larger 
than the thermal energy k^,T, where &b is the Boltzmann 
constant and T is the temperature, Vij(r) is approxi- 
mated as a quadratic form, 

%(r) ~ ~(r - R %j f = \[r- (R(v t ) + R{a 3 ))] 2 ■ (1) 
The condition of the applicability of this approximation 
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FIG. 2: LHS (Red circles) and RHS (Green triangles) of 
Eq. (4) vs 1/L. In the former, we put m A = —1 and tub = +1. 
In the latter, there is no restriction on tua and me. 



is discussed later. 

In this way, the Hamiltonian is given by 

N 2 L. N 

5 = Ey + 2Ete^^^) 2 + J} E' 7 - ( 2 ) 

i=l i=l 

Here <7i and pj are the coordinate and the momentum 
of ith. molecule. The symbol (i,j) means that a site j 
is one of the nearest neighbor sites of i. The last term 
represents effect of the ligand field. When D = 0, there 
are two ground states, i.e., <jj = I, Vi and <jj = —I, 
Mi. Throughout this Letter, according to the previous 
work [12], we fix the parameters as R(— 1) = 1, i?(l) = 
1.1, and k = 40. 

For convenience, we also label the molecules by the 
two-dimensional vectors r< — with Xi,yi G 

{1, 2, . . . , L}. The molecule at the xth column and yth 
row is labeled by (x,y), see Fig. 1(b). We distinguish 
from rii the former is a dynamical variable denoting the 
position of the ith molecule, but the latter is the label of 
the site denoting its positional relation on the lattice. 

Here we should mention the condition under which we 
can consider the equilibrium state of Eq. (2); the obser- 
vation time i bs should satisfy 

T cq < t ohs < ry oc exp[V /fc B T], (3) 

where r oq is the equilibration time of this model and ry 
is the activation time necessary to get over the potential 
depth Vq. If this condition were violated, fracture of the 
solid would occur and the approximation (1) is not neces- 
sarily valid; it becomes important to consider molecules 
whose positions are far away from the equilibrium posi- 
tions of the potential. When Vo » k B T, such slow things 
have not happened in the timescale i bs given above, and 
hence the model (2) is appropriate. 

The Hamiltonian (2) contains only the nearest- 
neighbor interactions and this model is extensive. Nev- 
ertheless, it is nonadditive, as we explain below. Let us 
consider the system in contact with a thermal bath at the 



temperature T. We virtually divide the system into two 
subsystems A and B; the subsystem A is composed of 
the molecules with xi < L/2 and the subsystem B with 
Xi > L/2. If the molecule i belongs to the subsystem A 
(B), we write i € A(B). We divide the Hamiltonian into 
three parts, H A b(X) = H A + Hb + XHj. Here Hx = 

E ie x(pf/2 + D^) + (*/2)E (w>eX (l« - *i\ - Rv) 2 > 
where X is A or B. The parameter A controls the strength 
of the interaction between subsystems and we choose it 
so that H A b(X = 1) is equal to Eq. (2). 

The concept of additivity is related to the probability 
that the magnetization of A is equal to Ma and that of B 
is equal to Mb . In this Letter, the system is said to be ad- 
ditive if P AB (T, M a , M b ) ~ P A (T, M A )P B (T, M B ). Here 
Pab{T, Ma, Mb) denotes the probability that the mag- 
netization of A at the temperature T is equal to Ma and 
that of B is equal to Mb in a composite system .Hab(I)- 
On the other hand, Px(T, M x ) (X is A or B) denotes the 
probability that the magnetization of X at the tempera- 
ture T is equal to Mx in a decoupled system Hx- 

According to statistical mechanics, such probabili- 
ties are related to the thermodynamic functions with 
some restriction [13]. That is, P A b(T,M a ,M b ) = 
exp{[F AB (T) - Fab(T,Ma,M b )]/£;bT}, where F AB (T) 
is the free energy of the composite system H A b(1) at 
the temperature T and F A b(T, Ma, Mb) is that un- 
der the constraint that the magnetizations of A and B 
are fixed to be Ma and Mb, respectively. Similarly, 
Px(T,M x ) = exp{[F x (T) - F X {T, M x )}/k B T}, where 
Fx{T) is the free energy associated with the Hamiltonian 
Hx and F X (T, M x ) is that under the constraint that its 
magnetization is fixed to be Mx- 

Thus the condition of additivity is rewritten as 

F AB (T, M a , M b ) - F A (T, Ma) - F B (T, M B ) 

~ F AB {T) - F A (T) - F B {T) (4) 

for any values of Ma and Mb- In particular, when the 
RHS of Eq. (4) is zero, the system is said to be exten- 
sive. According to thermodynamics, the LHS of Eq. (4) 
is equal to the amount of work done by the system in a 
quasi-static (reversible) isothermal process of changing A 
from 1 to 0. During the process, Ma and Mb should be 
individually conserved. 

In Fig. 2, we calculated the work per molecule, w — 
— J (dH (A) / OX) Xdt / L d , in such a quasi-static isothermal 
process. In the red circles, we put m A = M A /L d = — 1, 
rn B = M B /L d = +1, and T = 0.1. The value of w cor- 
responds to the LHS of Eq. (4). It is observed that w 
does not vanish as L increases. On the other hand, the 
green triangles are measured values of w without restric- 
tion on the magnetizations, which correspond to the RHS 
of Eq. (4). In this case w tends to zero as L increases. 
Therefore, it is concluded that the system is extensive 
but nonadditive. 

Because the system is not additive, it is expected that 
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FIG. 3: Equilibrium quantities. The red circles and the green triangles represent the results for the elastic spin model and 
for the effective Ising model, respectively, (a) Absolute magnetization densities, (|m|), are plotted against the temperature in 
the canonical ensemble, (b) Magnetic fields are plotted against magnetization densities at T = 0.26 in the restricted canonical 
ensemble, (c) Temperatures are plotted against energy densities in the microcanonical ensemble. In order to compare the two 
models, the contribution of the lattice vibration, 2T, is subtracted from the energy densities for the elastic spin model. 



the elastic spin model exhibits the peculiar properties ob- 
served in long-range interacting systems. We performed 
the Monte Carlo simulations to investigate equilibrium 
properties of the elastic spin model for (a) the canoni- 
cal ensemble, (b) the restricted canonical ensemble (the 
canonical ensemble with restriction on the value of the 
magnetization), and (c) the microcnaonical ensemble. 
Numerical results are shown by the red circles in Fig. 3. 
In the microcanonical ensemble, we subtracted 2T, which 
represents the contribution of the lattice vibration, from 
the total energy density. The parameters are set to be 
D = and g = 1 for (a) and (b), and D = 0.15, g = 20 
for (c).[19] The system size is L = 40. 

In the canonical ensemble, spin variables are changed 
according to the usual Metropolis algorithm and the po- 
sitions of molecules move according to the Hamilton dy- 
namics (by the leap-flog algorithm). We can see that 
the system undergoes a second order phase transition at 
T c ~ 0.35 in Fig. 3(a). The critical behavior belongs to 
the mean-field universality class [12] (see also Ref. [14]). 

The algorithm used in the restricted canonical ensem- 
ble is the same as in the canonical ensemble except that 
we prepare an excess degree of freedom referred to as the 
"demon" . The demon keeps the magnetization = ± 1 . 
Only if m&Oi < 0, the spin flip is accepted according to 
the Metropolis transition probability. After the flip, we 
change o~i and to — <7j and -mj, respectively. The ad- 
vantage of this method is that we can easily measure the 
magnetic field h = dF(T, M) jdM from the average value 
of ma by h = In prfe^j ■ In Fig- 3(b), we can clearly 
see the region where the susceptibility \ = dm / dh is neg- 
ative. The susceptibility is always positive in the canoni- 
cal ensemble and this discrepancy shows that the canon- 
ical ensemble is inequivalent to the restricted canonical 
ensemble. 

In the microcanonical ensemble, we used the leap- 
flog algorithm for the time evolution of {qi,Pi} and the 
Creutz algorithm [15] for the spin flip, in which an excess 



degree of freedom also called the "demon" , which has a 
positive energy, is prepared and the total energy of the 
system and the demon is conserved. The temperature of 
the system can be measured simply by T = (Ed). Fig- 
ure 3(c) clearly demonstrates that there is a region of 
the negative specific heat, which also shows the ensemble 
inequivalence. 

Although there is no direct interaction between a% and 
cjj in long distance, it will be a plausible consideration 
that an effective interaction arises between spin variables 
via lattice distortion. Nonadditivity of the system implies 
that the effective interaction is a sort of long-range inter- 
actions. In general, the effective Hamiltonian for {o~i} 
obtained by eliminating {qi,Pi} contains many-body in- 
teractions and will be very complicated, and hence it 
is hard to determine its explicit form. Here we assume 
that the effective interaction is a two-body interaction 
and try to estimate it from the numerical data of corre- 
lation functions. We put g = 1 and D — 0, and consider 
the canonical ensemble. The effective Hamiltonian is as- 
sumed to be written as H c g — —(1/2)5^ JijOiOj. We 
consider the high-temperature phase (T = 0.5). In this 
case, we can show that the relation between the correla- 
tion function C'ij = (o~i<jf) and the interaction potential 
is given by 



C(L)= I-J(L)/(k B T) 



(5) 



where the identity matrix is denoted by I. The rela- 
tion (5) is approximate one in general, but it becomes 
exact if the effective interaction is long-ranged, i.e., 



Jij(L) — L~ d (f>(L~ 1 r i j), 



(6) 



where r,-j = (xij , yij) = ri — Vj and d is the spatial dimen- 
sion (here d — 2). It is assumed that f^ <s <fi(x)d d x < 
+oo for an arbitrary fixed 6 > 0. It includes the power- 
law interaction <fi(x) ~ l/x a with a < d. The scaling 
L~ l Tij means that the interaction range is comparable 
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FIG. 4: (a) A global image of the interaction potential at 
T = 0.5 and L = 60. The site j is fixed at the center of the 
system (xj = yj = 30). In the yellow region Jij > 0.0003 
and in the black region Jij < —0.0003. (b) Graphs of the 
estimated effective interactions along the diagonal direction, 



for L = 40, 60, 



and 100. The transverse axis 



is the scaled distance Xij/L and the longitudinal axis is the 
scaled interaction matrix L 2 Jij. Only the points with Jy > 
are plotted. 



with the system size L. The scaling L~ d in Eq. (6) en- 
sures that the energy is proportional to the system size, 
which corresponds to so called Kac's prescription [3, 4]. 
The derivation of Eq. (5) will be reported elsewhere. 

In Fig. 4(a), the structure of the estimated interaction 
matrix Jij is depicted for L = 60. In the figure, the cen- 
tral site Tj = (L/2, L/2) is chosen as the site j. It is found 
that the interaction is highly anisotropic and a long-range 
ferromagnetic interaction emerges. Importantly, the spa- 
tial average of the interaction potential does not vanish; 
the long-range ferromagnetic interaction is not screened. 

The characteristic system size dependence of the esti- 
mated interaction is shown in Fig. 4(b), where the graphs 
of L d Jij(L) along the diagonal direction, = L — Xij, 
are depicted as a function of L~ 1 x i j for several values of 
L only in the region of Jij > 0. We can see that these 
graphs are collapsed well into a single curve (and this col- 
lapse occurs for any direction of ) . It means that the 
effective interaction is actually of the form of Eq. (6) and 
the use of Eq. (5) is indeed justified. Usually, Kac's pre- 
scription is regarded as a purely mathematical operation 
to extract nontrivial thermodynamic properties of long- 
range interacting systems [3]. However, in the present 
model, such a mathematical operation is not necessary; 
the scaling of Eq. (6) naturally emerges in the effective 
interaction. This is a remarkable characteristic of this 
model. 

In Fig. 3, we compared some equilibrium quantities of 
the elastic model with those of the effective Ising model, 
-ffoff = — ^2ij Jij&icrj — D^2,f Pi- In all the simulations 
presented in Fig. 3, we used the effective interaction J^ 
estimated at T — 0.5 in the canonical ensemble. The nu- 
merical results of the effective Ising model are plotted by 
green triangles. In all the three graphs, the results of the 



original elastic model excellently agree to those of the ef- 
fective Ising model. Although the effective interaction is 
estimated by the data at the single point in the canonical 
ensemble, the elastic model is indistinguishable from the 
effective Ising model for all the values of parameters and 
for all the statistical ensembles. 

To summarize, it was found that the Hamiltonian of 
the elastic spin model is local and extensive but nonad- 
ditive. The elastic spin model exhibits peculiarities like 
the negative specific heat and the negative susceptibility 
according to the ensemble. In addition, we numerically 
showed that as far as equilibrium states of spin degrees of 
freedom are concerned, the elastic spin model is indistin- 
guishable from the effective Ising model with a long-range 
interaction. It implies that statistical mechanics of long- 
range interacting systems is relevant for a certain kind 
of short-range interacting systems. The limitation of the 
present study is that several timescales should be sepa- 
rated as Eq. (3), which requires the large potential depth 
Vq compared to fceT; otherwise considering an equilib- 
rium state of the Hamiltonian (2) would not be justified. 
In the original system without the approximation (1), 
the system may restore additivity after a very long time 
lobs 3> Ty by occurring spontaneous fracture of the crys- 
tal. Finally, it is pointed out that the present study will 
give some insight into experimental attempts to realize 
a long-range interacting system and observe its peculiar 
properties in laboratory [16, 17]. Such attempts are not 
satisfactory yet, but we hope that this work stimulates 
different experimental approaches toward it. 
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